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Executive  Summary 


The  long-term  goal  of  this  research  is  to  develop  algorithms  to  simulate  high  Reynolds 
number  turbulent  flow  in  complicated  geometries  using  embedded  boundary  grids.  The 
main  stumbling  block  is  the  lack  of  an  affordable  refinement  strategy  that  is  common¬ 
place  in  structured  grids  -  namely,  the  use  of  highly  refined,  high  aspect  ratio  cells  in  the 
boundary  layer.  Over  the  grant  period,  we  developed  a  two-dimensional  viscous  steady- 
state  flow  solver  for  both  laminar  and  turbulent  flow.  We  developed  a  way  to  incorporate 
wall  functions  in  non-body  fitted  grids  by  having  each  cut  cell  spawn  a  “linelet”  from 
the  wall  through  the  cut  cell  centroid,  into  the  Cartesian  grid.  We  developed  a  fully 
conservative  method  for  coupling  the  wall  function  to  the  flow  on  the  Cartesian  grid. 
We  extended  the  wall  function  into  a  more  general  wall  model  that  solves  a  two  point 
boundary  value  problem  on  the  linelets.  The  bvp  includes  more  terms  than  is  typically 
used  in  the  diffusion  model,  which  is  what  standard  wall  functions  are  based  on.  The 
improved  wall  model  on  the  linelets  allow  them  to  go  further  into  the  boundary  layer, 
reducing  the  Cartesian  resolution  requirements  and  allowing  coarser  near- wall  Cartesian 
cells.  This  work  was  in  collaboration  with  Michael  Aftosmis,  at  NASA  Ames  Research 
Center. 


Research  During  Grant  Period 


In  previous  work  during  the  last  grant  period  we  developed  Navier  Stokes  solvers  for 
laminar  and  viscous  flows  in  two  dimensions.  In  integral  form  the  two  dimensional 
Reynolds- averaged  Navier  Stokes  equations  are: 


d 

dt 


U  dA+  d>  (Fi  +  Gj)  •  hdS  =  (f>  ( Fvi  +  Gvj ) 


aft 


n< 


aft 
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The  Prandtl  number  is  given  by  Pr  =  pcp/k ,  and  p  is  computed  using  Sutherland’s 


law  [8].  We  take  Pr  =  .72.  The  shear  stresses  are 


txx  =  -  (p  +  pt){2ux  -  Vy) 

Txy  =  (p  +  Pt)(uy  +  VX)  =  TyX  (2) 

2 

Tyy  =  g  (p  +  pt){2vy  -  ux) 

and  turbulent  heat  flux  Vq  =  +  ^^)VT. 

For  laminar  flow,  we  compared  several  stencils  for  use  at  cut  cells,  found  two  with 
roughly  equivalent  accuracy,  and  implemented  the  one  that  fit  better  in  Cart3D’s  finite 
volume  infrastructure.  To  be  able  to  compute  skin  friction,  which  needs  du/dy  at  the 
wall,  we  used  a  quadratic  in  the  normal  direction  in  the  cut  cells,  instead  of  the  linear 
reconstruction  used  in  the  rest  of  the  mesh. 

For  turbulent  flow,  we  used  the  RANS  equations  and  the  Spalart-Allmaras  (SA) 
turbulence  model.  The  SA  turbulence  model  defines  the  pt  =  P^t,  the  turbulent  viscosity 
in  Eq.(2): 
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Standard  forms  for  these  terms  are 

1  X>0 
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and  S  =  ,  where  ft  is  the  magnitude  of  vorticity.  For  the  destruction  term  we  have 


Destruction  = 
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with  i/t  =  vfv i  and  the  usual  definitions  of  fw ,  g ,  r,  x,  and  the  other  constants [13].  We 
assume  the  flow  is  fully  turbulent  and  do  not  include  transition  terms. 

We  developed  a  somewhat  modified  “negative”  model  to  handle  the  strong  transient 
behavior  on  coarse  grids  using  multigrid.  ( Negative  refers  to  how  to  handle  the  non¬ 
physical  negative  viscosity;  at  convergence  when  viscosity  is  back  to  being  positive  the 
equations  are  the  same).  Many  researchers  use  a  first  order  finite  volume  scheme  to 
compute  in  (3)  above,  to  help  keep  the  viscosity  positive.  We  found  it  necessary  to 
use  a  second  order  method  for  i> ,  since  on  non-flow-aligned  grids  this  was  simply  too 
diffusive. 

We  developed  a  hierarchy  of  wall  models  that  improve  on  the  wall  function,  allowing 
us  to  use  coarser  background  grids.  The  previous  research  period,  where  we  developed 
the  “quadratic”  wall  model  used  for  laminar  flow,  and  the  viscous  stencil  used  in  the 
highly  irregular  cut  cells,  has  previously  been  reported,  and  is  already  in  the  literature 
[2].  In  this  research  period  we  developed  a  wall  function  and  wall  model  approach, 
summarized  below. 


Analytic  Wall  Functions  for  Embedded  Boundary  Grids 

For  RANS  simulations,  the  near- wall  flow  is  governed  by  viscous  stresses,  as  opposed 
to  the  outer  flow  where  inertial  forces  dominate.  This  is  the  basis  for  the  use  of  wall 
functions,  which  have  been  receiving  renewed  interest  for  both  body-fitted  [6,  12,  9,  11] 
and  non-body-fitted  immersed  boundary  grids  [4,  7].  By  looking  at  the  asymptotics  of 
the  RANS  equations,  keeping  the  largest  terms,  and  assuming  for  simplicity  there  is  no 
pressure  gradient,  the  so-called  “diffusion  model”  can  be  derived, 

SM)- 

In  boundary  layer  (“plus”)  coordinates  and  after  integrating  once  this  becomes 

(1+^  =  1  <8> 


where 


z/t+  =  vt/v,  U+  =  U/uT ,  and  y+  =  yuT/v . 


There  are  two  regimes,  ut+  <  1,  which  gives  the  linear  behavior  U+  =  y+  in  the 
viscous  sublayer,  and  vf  1,  which  gives  the  logarithmic  layer  U+  =  (1  /«)  log (y+)  +  B, 
with  k  =  .41,  when  the  turbulent  viscosity  is  approximated  using  The 

constant  B  =  5.0333.  Spalding’s  wall  function,  given  by 

y+(u+ )  =  u+  +  e~KB  -  1  -  ku+  -  ^{k,u+)2  -  i(/s«+)3^  (10) 

is  a  functional  fit  to  these  two  regimes. 
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Figure  1:  Comparison  of  SA  wall  model  with  Spalding’s  formula.  The  largest  difference  is 
in  the  transition  region. 


Recently  there  is  a  new  wall  function,  developed  by  Allmaras  [2,  1],  which  is  more 
complicated  to  derive  but  more  compatible  with  the  limiting  form  of  the  Spalart- Allmaras 
model  used  in  our  RANS  simulations,  see  Fig.  1.  In  final  form  it  is 

u+{y+)  =  B  +  c\  log((y+  +  ai)2  +  b\)  -  c2  log((y+  +  a2)2  +  b\) 

—  C3ArcTan(y+  +  ai,  bi)  —  C4ArcTan(y+  +  a2,b2), 

where  we  omit  the  values  for  the  many  constants.  Equation  (11)  (henceforth  called 
the  SA  wall  model)  has  the  advantage  over  Spalding’s  formula  of  matching  the  pro¬ 
files  actually  computed  in  the  flow  field  by  the  Spalart- Allmaras  turbulence  model.  In 
particular  the  profiles  are  a  better  match  through  the  transition  region,  where  the  two 
formulas  differ  the  most.  This  is  evident  in  Fig.  1.  In  this  figure,  Spalding’s  formula  was 
evaluated  using  B  =  5.033  in  (10)  to  match  the  asymptotic  behavior  of  the  SA  model, 
which  produces  this  value  for  the  constant  shift.  Furthermore,  (11)  is  more  computa¬ 
tionally  efficient  since  given  y+ ,  the  value  u+  is  explicitly  evaluated  rather  than  needing 
a  nonlinear  iteration  as  (10)  does.1 

The  wall  function  is  coupled  to  the  underlying  Cartesian  grid  through  its  endpoints. 
This  is  illustrated  schematically  in  Fig.  2.  At  the  wall  it  is  anchored  with  a  zero  velocity 
boundary  condition  at  the  wall  in  each  cut  cell  C.  At  the  other  end,  there  is  a  fixed 
point  F  (to  use  the  same  terminology  as  [4])  located  a  distance  hp  from  the  boundary 
in  the  normal  direction,  which  anchors  it  in  the  flow.  (The  use  of  a  fixed  distance  away 
from  the  wall,  instead  of  the  cut  cell  centroid  for  example  is  needed  to  compensate 
for  the  wildly  irregular  cut  cells  in  the  mesh.)  The  solution  at  point  F  is  obtained  by 
interpolation  from  the  Cartesian  mesh.  The  velocity  is  then  rotated  into  a  normal  and 
tangential  frame  of  reference  to  obtain  qtang ,  the  tangential  velocity,  using  the  normal 
vector  defined  by  the  portion  of  the  boundary  in  cut  cell  C.  A  Newton  iteration  is  used 

1  However  Kalitzin  et  al.  [9]  have  a  nice  approach  that  avoids  iteration  using  pre-computed  tables  of 
inverses. 


4 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


Figure  2:  Illustration  of  construction  for  wall  model  coupling  via  constant  height  “forcing 
pts”  through  the  centroids  of  the  cut  cells. 


to  find  the  friction  velocity  uT  (see  eq.  (9)  )  so  that  the  points  lie  on  the  wall  function. 
Very  few  iterations  are  needed  since  the  uT  from  the  previous  timestep  provides  a  great 
starting  guess  for  the  next  iteration. 

Gradients  of  the  model  provide  the  viscous  flux  du/dy\y=o  needed  at  the  wall  by 
the  finite  volume  scheme.  The  model  also  gives  values  of  the  tangential  velocity  needed 
to  compute  the  difference  at  the  cut  faces  (after  rotating  back  to  Cartesian  velocity 
components  u  and  u).  Finally,  we  improve  the  gradient  originally  computed  in  a  cut  cell 
with  the  gradient  of  the  solution  to  the  wall  model,  evaluated  at  the  cut  cell  centroid. 
(This  was  also  done  with  the  quadratic  model  for  laminar  flow).  This  provides  a  fully 
conservative  coupling  of  the  wall  function  and  mesh,  and  makes  maximal  use  of  the  wall 
function. 

Figure  3  shows  computational  results  from  turbulent  flow  over  a  bump  in  a  channel, 
taken  from  the  NASA  Langley  Turbulence  Modeling  Resource  [14]  so  we  can  compare 
with  their  results  using  CFL3D. 

BVP  Wall  Models 

There  are  several  drawbacks  to  wall  functions.  For  flows  with  separation,  the  skin  fric¬ 
tion,  and  wr,  go  to  zero  and  the  flow  reverses.  We  will  not  be  able  to  compute  separation 
using  a  model  in  plus  coordinates  since  that  entails  dividing  by  zero.  Furthermore,  the 
diffusion  model  assumes  the  convective  terms  are  small,  so  cannot  be  used  far  from  the 
wall  where  this  is  no  longer  true.  It  also  assumes  a  mixing  length  model  [10]  for  the 
turbulent  viscosity,  which  is  not  true  away  from  the  wall  after  the  initial  linear  growth. 

We  have  developed  a  hierarchy  of  two  point  boundary  value  problems  (bvp)  that 
include  increasingly  more  “physics” ,  whose  solution  replaces  the  wall  function.  The  first 
step  was  to  actually  replicate  the  wall  function  as  a  two  point  bvp.  The  diffusion  model 
wall  function  (7)  is  already  a  bvp.  By  using  the  mixing  length  model  for  /i*  it  can  be 
easily  solved  and  tested  in  our  framework.  We  use  the  Shampine  BVP  solver  for  this, 
modified  to  run  in  parallel  in  a  thread-safe  way,  save  the  solution  between  iterations, 
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skin  friction  <SA) 


Figure  3:  Skin  friction  using  SA  wall  model  on  turbulent  bump  in  channel  problem.  The 
wall  model  stays  the  same  as  the  grid  is  refined.  Without  the  wall  model,  many  more  mesh 
refinements  are  needed  beyond  our  finest  level  grid. 


etc.  The  solver  takes  a  systems  of  first  order  equations,  in  this  first  case  it  is  two 
equations,  and  u( 0)  and  u(F)  for  the  boundary  conditions.  As  before,  u(F)  is  found  by 
interpolation  from  the  Cartesian  grid.  This  eliminates  the  problem  of  uT  — v  0,  since  this 
works  in  physical  coordinates  and  not  plus  coordinates.  We  then  evaluate  the  solution 
du/dy\y=o  to  get  uT,  and  compute  the  wall  boundary  conditions  on  the  Cartesian  side. 
The  mixing  length  model  we  use  is 

vt  =  Xv  =  Ky+v  and  y+  =  yuTReL/M (12) 
so  that  the  turbulent  viscosity  is  proportional  to  the  distance  from  the  wall. 

p 

mixing  length  model  ....  mixing  length  model 


Figure  4:  Results  showing  very  different  behavior  of  eddy  viscosity  profiles  at  different 
stations  along  wall,  taken  from  very  fine  reference  solution.  AGARD  case  10  RAE  2822 
airfoil  simulation  [5]. 


However  as  Fig.  4  illustrates,  this  is  only  the  first  part  of  the  turbulent  viscosity 
computed  by  the  RANS  equations,  valid  only  very  close  to  the  wall. 


6 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


Figure  5:  Results  of  solving  the  bvp  using  a  mixing  length  model  and  an  SA  model  for 
viscosity.  The  curvature  in  the  eddy  viscosity  is  well  captured,  and  the  tangnetial  velocity 
changes  somewhat  as  well. 


The  next  step  up  in  the  hierarchy  is  to  replace  the  mixing  length  model,  and  use 
as  much  of  the  Spalart-Allmaras  turbulence  model  for  vsa  as  is  locally  available.  For 
example,  all  the  normal  derivatives  that  are  already  part  of  the  bvp  are  included.  We 
do  not  compute  the  terms  involving  streamwise  derivatives.2.  This  gives  the  model 


s' 

d_ 

dy 


du 


dv 

dy 


dp 

dx 

=  rest  of  diffusion  terms  +  Prod  +  Dest 


(13) 

(14) 


When  written  as  a  system  of  first  order  equations  this  has  dimension  4,  and  the  new 
boundary  conditions  are  vt( 0)  =  0  at  the  wall,  and  Vt (-F) ,  again  specified  via  interpolation 
from  the  Cartesian  grid.  We  initialize  the  solution  using  the  previous  iterate,  or  initially 
using  the  SA  wall  function  and  a  mixing  length  model  using  the  SA  turbulence  model. 
This  helps  the  BVP  solver  converge  quickly. 

Using  (most  of)  the  SA  turbulence  equations  allows  much  more  variation  in  turbu¬ 
lent  viscosity  to  be  computed,  see  e.g.  the  profiles  in  the  bottom  row  of  Fig.  4.  The 
streamwise  velocity  and  especially  the  turbulent  viscosity  profiles  vary  greatly,  depend¬ 
ing  on  the  location  along  the  airfoil  and  the  boundary  layer  thickness  there.  The  four 
equation  bvp  model  is  able  to  compute  the  eddy  viscosity  even  in  this  very  nonlinear 
regime.  Figure  5  shows  a  computation  with  our  model  compared  to  the  mixing  length 
model.  The  curvature  in  the  eddy  viscosity  profile  is  captured  very  well. 

We  experimented  for  a  long  time  with  computing  streamwise  derivatives  of  terms  in 

2  We  use  standard  terminology  and  refer  generically  to  the  u  velocity  as  the  streamwise  velocity,  even 
though  with  cut  cells  it  is  more  complicated,  similarly  for  the  v  velocity,  instead  of  saying  normal  and 
tangential  velocities 
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these  models,  but  found  this  to  be  too  noisy  to  work.  This  involves  differentiating  from 
stations  at  the  same  normal  location  on  each  linelet.  However  the  bvps  are  only  solved 
to  a  tolerance,  and  differentiating  them  amplifies  the  noise.  There  is  also  noise  coming 
from  the  changing  location  of  point  F  relative  to  the  Cartesian  mesh,  for  example  in  the 
interpolation  errors.  Finally,  there  is  the  problem  with  the  equation  for  the  v  velocity  - 
the  velocity  in  the  normal  direction.  Asymptotically,  this  is  really  an  equation  for  the 
pressure,  and  since  it  is  orders  of  magnitude  smaller  than  the  streamwise  velocity  it  is 
very  hard  to  compute. 

Instead,  our  highest  fidelity  model  adds  some  convective  terms  taken  from  the  back¬ 
ground  Cartesian  grid  to  the  equation  for  the  tangential  velocity.  Those  terms  are  zero 
at  the  wall,  but  as  you  go  further  into  the  boundary  layer  and  approach  the  inviscid 
region,  they  are  needed  to  balance  the  pressure  gradient.  We  get  the  convective  terms 
at  point  F  again  by  interpolation  from  the  Cartesian  grid.  However  we  have  to  be  able 
to  evaluate  them  at  any  point  y  G  [0,  F]  to  solve  the  bvp.  For  this  we  use  a  shut  off 
function  that  goes  from  0  at  the  wall,  to  1  at  point  F.  Altogether,  our  enhanced  wall 
model  together  with  the  equation  for  the  turbulent  viscosity,  is 


d_ 

dy 


((v  +  ih) 


d_ 

dy 


\v  +  v) 


=  1 


.  du(F)  .  du(F) 


dx  dy  \ 

=  rest  of  diffusion  terms  +  Prod  +  Dest 


(15) 


Our  model  for  ip  is 

1%)  =  SA(y)/SA(F), 

which  satisfies  ip(0)  =  0  and  ip(F)  =  1,  where  SA  is  the  Spalart-Allmaras  wall  function. 

We  are  submitting  a  paper  with  the  details  on  this  wall  model  to  the  AIAA  in  the 
coming  month.  Here  we  just  show  one  result  using  this  approach.  We  use  an  example 
of  flow  around  a  NACA  0012  from  the  Turbulence  Modeling  Resource  website.  There 
is  a  strong  pressure  gradient  around  the  leading  edge.  Figure  6  shows  skin  friction  for 
both  the  SA  wall  function  and  the  bvp  wall  model.  The  SA  wall  function  is  not  properly 
centered  with  respect  to  the  converged  CFL3D  solution  on  either  mesh,  whereas  the  bvp 
model  is  aligned  on  even  the  coarser  of  the  two  mesh  results.  Other  researchers  have 
tried  without  success  to  include  only  a  pressure  gradient  on  the  right-hand-side  of  the 
diffusion  model  [3].  We  believe  that  is  because  the  convective  terms  are  of  the  same 
order  as  the  pressure  gradient,  and  all  terms  are  needed  to  find  the  right  balance.  Note 
that  the  coarser  of  the  two  meshes  in  Fig.  6  has  y+  «  500  at  the  leading  edge! 
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Figure  6:  Comparison  of  SA  wall  function  and  bvp  wall  model  for  the  NACA0012  case.  This 
example  convincingly  shows  how  including  the  pressure  gradient  term  centers  the  profile. 
The  grid  labeled  'Fine'  (as  opposed  to  Veryfine)  has  y+  ~  500  at  the  leading  edge. 
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